﻿#include "opencv2/opencv.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"

#include <iostream>
#include <fstream>
#include <numeric>

using namespace cv;
using namespace std;

//4. 使用米粒图像，分割得到各米粒，首先计算各区域(米粒)的面积、长度等信息，进一步计算面积、长度的均值及方差，分析落在3sigma范围内米粒的数量。 

int main(void)
{
	char *fn = (char *)"E:\\AI\\AIsoftware\\opencv\\sources\\samples\\data\\rice.png";
	Mat image = imread(fn);
	//imshow("原图",image);

	Mat gray, bw; //二值化后的图像

	//灰度图转换
	cvtColor(image, gray, COLOR_BGR2GRAY);
	//大津算法阈值化
	threshold(gray, bw, 0, 0xff, CV_THRESH_OTSU);
	imshow("阈值化", bw);

	//形态学处理，去噪
	Mat element = getStructuringElement(MORPH_CROSS, Size(3, 3));  //返回指定形状和尺寸的结构元素,交叉形：MORPH_CORSS;
	//morphologyEx(bw, bw, CV_MOP_CLOSE, element);  //形态学变换函数，闭运算，先膨胀后腐蚀
	morphologyEx(bw, bw, CV_MOP_OPEN, element);//开运算
	
	//图像分割
	Mat seg = bw.clone();  //复制一份
	vector<vector<Point>> cnts;	
	findContours(seg, cnts, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);  //提取轮廓

	/*
	//绘制轮廓
	vector<vector<Point>> hierachy;
	Mat draw = Mat::zeros(seg.size(), CV_8UC1);
	for (size_t i = 0; i < cnts.size(); i++)
	{
		drawContours(draw, cnts, i, (200, 200, 200), 1, 8, hierachy, 0, Point(0, 0));
	}
	imshow("轮廓", draw);
	*/

	//筛选
	int count = 0;      //米粒总数
	set<double> area_set;//string集合
	set<double> length_set;//string集合

	//double a[300];
	//double le[300];
	//double arv_area=0.0;
	//double arv_length=0.0;
	string strCount;

	for (int i = cnts.size()-1; i >0; i--)
	{
		vector<Point> c = cnts[i];
		double area = contourArea(c);  //计算整个或部分轮廓的面积
		double length = arcLength(cnts[i], true); //长度
		if (area < 10)  //滤除面积小于10的分割结果：可能是噪声
			continue;

		//统计数量
		//arv_area += area;
		//arv_length += length;
		area_set.insert(area);
		length_set.insert(length);
		count += 1;
		cout << "blob" << i << " 面积:" << area<<" 长度："<<length<<endl;   //打印出每个米粒的面积、长度

		Rect rect = boundingRect(c);     //提取矩形坐标，boundingRect()计算轮廓的垂直边界最小矩形
		//cout << " x:" << rect.x<< " y:" << rect.y << endl;

		//在原始图像上画出包围矩形，并给每个矩形标号
		rectangle(image, rect, Scalar(0, 0, 255), 1);
		stringstream ss;
		ss << count;
		strCount = ss.str();
		putText(image, strCount, Point(rect.x, rect.y), CV_FONT_HERSHEY_PLAIN, 0.5, Scalar(0, 255, 0)); //在图像上绘制文字
	}
	imshow("imageshow", image);
	//waitKey(0);

	//米粒面积的均值和方差
	double sum1 = accumulate(std::begin(area_set), std::end(area_set), 0.0);
	double averageArea = sum1 / area_set.size(); //均值

	double accum1 = 0.0;
	std::for_each(begin(area_set), end(area_set), [&](const double d) {
		accum1 += (d - averageArea)*(d - averageArea);
	});
	double stdevArea = sqrt(accum1/(area_set.size() - 1)); //方差

	//米粒周长的均值和方差
	double sum2 = accumulate(std::begin(length_set), std::end(length_set), 0.0);
	double averageLegth = sum2/area_set.size(); //均值

	double accum2 = 0.0;
	std::for_each(begin(length_set), end(length_set), [&](const double d){
		accum2 += (d - averageLegth)*(d - averageLegth);
	});
	double stdevLength = sqrt(accum2/(length_set.size() - 1)); //方差

	cout << "米粒总数:" << count << endl;	
	cout << "平均面积：" << averageArea << " 平均长度：" << averageLegth << endl;
	cout << "面积方差:" << stdevArea << " 长度方差:" << stdevLength << endl;

	//落在3sigma范围内米粒的数量
	double sigmaArea = sqrt(stdevArea);
	int areaIn3sigma = 0;	
	std::for_each(begin(area_set), end(area_set), [&](const int d){
		//areaIn3sigma = d - averageArea;
		if ((d - averageArea)<=3* sigmaArea)
			areaIn3sigma += 1;
	});

	double sigmaLength = sqrt(stdevLength);
	int lengthIn3sigma = 0;
	std::for_each(begin(length_set), end(length_set), [&](const int d) {
		//lengthIn3sigma = d - averageLegth;
		if ((d - averageLegth) <= 3 * sigmaLength)
			lengthIn3sigma += 1;
	});
	cout << "落在3sigma范围内米粒的数量（面积）:" << areaIn3sigma << " 落在3sigma范围内米粒的数量（长度）:" << lengthIn3sigma << endl;

	waitKey(0);
}